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Abstract 

A computational scheme for solving 2D Laplace boundary-value problems using rational 
functions as the basis functions is described. The scheme belongs to the class of desin- 
gularized methods, for which the location of singularities and testing points is a major 
issue that is addressed by the proposed scheme, in the context of the 2D Laplace equa- 
tion. Well-established rational-function fitting techniques are used to set the poles, while 
residues are determined by enforcing the boundary conditions in the least-squares sense 
at the nodes of rational Gauss-Chebyshev quadrature rules. Numerical results show that 
errors approaching the machine epsilon can be obtained for sharp and almost sharp cor- 
ners, nearly-touching boundaries, and almost-singular boundary data. We show various 
examples of these cases in which the method yields compact solutions, requiring fewer 
basis functions than the Nystrom method, for the same accuracy. A scheme for solving 
fairly large-scale problems is also presented. 



1. Introduction 

Laplace boundary-value problems (BVPs) are prevalent in the applied sciences, and 
methods for their solution have been the subject of study for decades. Among these, 
prominent approaches are finite-difference and finite-element methods, integral-equation 
methods, and, in 2D, conformal-mapping techniques. While the Laplace equation is the 
simplest elliptic partial differential equation (PDE), improving methods for its solution 
is nevertheless of continued interest. At the focus of recent research are problems with 
small radii of curvature, nearly-touching boundaries, and singular, or almost-singular 
boundary data ^nQ]- 

The technique described in this paper belongs to a class of methods known as desin- 
gularized methods, which are used for the solution of linear BVPs. What characterizes 
these methods is that the basis functions are elementary solutions of the homogenous 
PDE, and hence the singularities of the basis functions lie outside the EVP domain. 
The simplest example of a desingularized method is the method of images, which is only 
applicable in a handful of special cases for which the exact locations and amplitudes of 
the image sources are known analytically. For more general problems, an approximate 
solution can be obtained by setting the locations of the sources heuristically and solv- 
ing for the amplitudes so that the boundary conditions are satisfied at a set of testing 
points. Various types of image sources and heuristics have been proposed over the years. 
Vekua ;2] and Kupradze Q are usually credited as the pioneers of these methods, and 
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other early works include those of Fox, Henrici, and Moler ^ , Bergman [l^] , and Eisen- 
stat [llj. Nowadays, the methods in this category go by various names, such as method 
of fundamental solutions 12|, |l3| , method of auxiliary sources [l5| . generalized mul- 
tipole tech niq ue [l6| . method of fictitious sources Ch. 1.5.6], and source-model 

technique 19l-l21|. 



Desingularized methods share many of the properties of integral-equation methods, 
and indeed can be viewed as such when the sources are distributed on a curve lying 
outside of the BVP domain. The basis functions used for the density functions are then 
Dirac delta functions. As the potentials due to these basis functions are smooth on the 
boundaries, it is customary to enforce the boundary conditions at a set of testing points. 
This obviates laborious panel integrals that slow down conventional Galerkin schemes 
for boundary integral equations, making the cost of forming the system matrix similar 
to that of a Nystrom method. Also, the representation of the solution with elementary 
sources is very convenient for post-processing tasks such as calculating capacitances, 
charges, forces, etc., and the potential itself at arbitrary points. 

Owing to the above properties, desingularized methods are among the simplest meth- 
ods for solving linear EVPs. They can also be very accurate, provided the sources and 
testing points are located "adequately" . To substantiate this claim, we show a couple of 
results obtained with the method proposed in this paper. In Fig. [TJi, 30 dipoles are used 
to solve the Laplace equation outside two closely spaced curves. The solution is required 




Figure 1: Solutions obtained with a small number of basis functions, for (a) exterior and 
(b) interior Laplace Dirichlet problems. The number of dipoles is 30 in (a) and 60 in (b). 
Darker dots correspond to larger dipole amplitudes. 



to be —1 on the lower curve and +1 on the upper curve, and this boundary condition is 
satisfied with a maximum error of A£',„ax — 1-7 x 10^^. In Fig. [TJd, an interior Dirich- 
let problem with two perfectly sharp reentrant corners is solved using 60 dipoles. The 
boundary condition is due to a monopole located inside the polygon and the maximum 
error is A_Emax — 2.5 x 10~^. The Nystrom method, which is spectrally accurate when 



the boundaries are smooth [24I, is the standard against which the proposed method is 



compared with in this paper. Such a comparison, shown in Fig. [^J reveals that when 

2 



the boundary curves nearly touch, a desingularized method can indeed be much more 
accurate than the Nystrom method. 




Figure 2: For the problem shown in Fig. [T^, maximum errors on the order of 10~^^ 
are obtained with no more than 70 dipoles. Using the Nystrom method with the same 
number of basis functions yields errors of about 5%. 



The adequate location of the sources and testing points is a major unresolved issue 
with desingularized methods, and much work has been devoted to its analysis ( (23l - [27j | 
are just a few references out of a large body of work). Source location affects both the 
rate of error decay and the effect of round-off errors which may limit the highest accu- 
racy achievable. Also, the effective application of the fast multipole method (FMM) [2^ 
appears to depend strongly on the source location (29j . The problem with applying the 
FMM to a desingularized method is that the linear systems to be solved are generally 
ill-conditioned. In [2^, the authors proposed a block-diagonal preconditioner that accel- 
erates convergence, provided the distance between the sources curve and the boundary 
curve is less than a problem dependent critical distance. 

Source location schemes that guarantee exponential convergence were developed in 
and [2^, [33, for Laplace and Helmholtz Dirichlet BVPs in simply-connected analytic 
domains, respectively. In these works, sources are distributed uniformly (in so-called 
conformal angle) on curves conformal with the boundary, and the resulting error in the 
solution is shown to decrease exponentially, in exact arithmetic. But placing sources on 
a conformal curve may be restrictive, and any a-priori source restriction requires a-priori 
knowledge of singularity locations. 

A common approach has been to formulate the source-location problem as an opti- 
mization problem. If the number of sources is set in advance, finding the best source 
locations is a nonlinear least-squares problem, and various attempts have been made 
to use nonlinear optimization techniques to determine the source locations. Unfortu- 
nately, this optimization problem is not convex and has many local minima. Therefore, 
descent-type methods such as Gauss-Newton and Levenberg-Marquardt are sensitive to 
the initial source locations [12', '31' J [HI, while more robust methods that have been tried, 
such as simulated annealing [33] and genetic algorithms 3J], may take a very long time 
to converge. Moreover, if the positions of the sources are adapted, the positions of the 

3 



testing points must be adapted as well. If this is not done, and the sources approach 
the boundary, the error at the testing points will not be representative of the error on 
the boundary. Alternatively, if the cost function of the optimization is defined as an 
integral on the squared error (as in [s^ . for example), the quadrature for this integral 
must be adapted along with the source positions. This issue of adapting the testing 
points has been largely ignored in the literature; in this paper it is addressed by use 
of rational Gauss-Chebyshev quad rature rules. Perhaps the most practical approaches 
use geometry-based heuristics |l6l [25l . [ssl [36| . While these can be more than adequate 
in many engineering applications, they fail to take full advantage of the flexibility in 
locating the sources and their accuracy may be limited. 

In this paper we describe a new method for determining appropriate locations for the 
sources and testing points. The method is tailored to the 2D Laplace equation, for which 
the solution can be represented as the real part of a complex potential which we ap- 
proximate by the sum of a rational function and (for multiply-connected domains) a few 
logarithmic terms. We use well-established techniques for finding rational-function ap- 
proximations to the boundary data, and the poles of these rational functions correspond 
then to the positions of dipole sources. 

The major difficulty in applying this idea is that, in a Dirichlet problem, the bound- 
ary data corresponds to the real part of the complex potential on the boundary. The 
imaginary part is determined, up to an additive constant, by the Dirichlet-to-Neumann 
map. Of course, applying this map is tantamount to solving the Dirichlet problem. Sim- 
ilarly, in a Neumann problem, only the imaginary part of the complex potential on the 
boundary is known. The solution we propose for this problem is to estimate the missing 
information (say, the imaginary part in a Dirichlet problem) by solving the BVP with a 
given set of poles and then using this information to improve the pole locations. Assum- 
ing the pole locations are indeed improved, they can be used to find a better estimate of 
the imaginary part, which in turn is used to improve their locations further. We proceed 
in this way until the pole positions appear to have converged, the error drops below a 
preset tolerance, or a maximum number of iterations is reached. 

Clearly, the process just described is not guaranteed to converge to the solution 
that would be obtained were the imaginary part known in advance. It is possible that 
the estimated imaginary part will be inaccurate, and the poles found based on this 
estimate will not yield any improvement on the next iteration. When this happens, the 
improvement in accuracy stalls. It is then usually, though not always, possible to increase 
the accuracy further by increasing the number of poles. In numerical experiments, we 
have found that despite this difficulty, errors approaching the machine epsilon can be 
obtained for various challenging configurations. 

2. Dirichlet problem, simply-connected domain 

We begin by considering the simplest case, in which we seek a real valued potential, 
U{x,y), which obeys the Laplace equation. More precisely, 

V^U{x, y) = 0, for all [x, y) G Q, (1) 

where O is a simply connected domain in bounded by a curve F, and U must also 
obey a Dirichlet boundary condition, 

C^lo„r = /- (2) 
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The method we propose for solving this problem is to approximate the solution as a 
weighted combination of dipole potentials, and to use rational function fitting algorithms 
to determine the dipole locations. 



2.1. Dipole representation 

As is well-known, the potential U can be sought as the real part of a complex function 
W{z), where z = x + iy and W(z) is holomorphic and single valued in il (it is single 
valued because Q is simply connected). The classical single- and double- layer potentials, 



WsinM^^) = / <J{Z') \0g{z' - Z)dz\ (3) 

and 



r 



W'doublc(^) - / 4^dz', (4) 

Z7^^ Jy' z' — z 

where a{z) and ^{z) are referred to as the monopole and dipole densities, respectively, 



have been used extensively to formulate integral equations [22|, |37|, |33| . For example, 
a second-kind Fredholm equation is readily obtained by combining the double-layer po- 
tential, with the Dirichlet boundary condition, ([2]), or by combining the single- layer 
potential, ([3]), with the Neumann boundary condition. 

In desingularized methods, the most common approach is to use the single-layer 
potential, 

Wv'iz) = I a{z')\og{z' -z)dz', (5) 



where F' encloses, but does not intersect, F. Enforcing the boundary condition on F 
leads to an ill-conditioned Fredholm equation of the first kind. Nevertheless, as discussed 
in [2^, if W{z) can be continued analytically throughout the region enclosed between 
F and F', highly accurate solutions can be obtained, despite this ill-conditioning. The 
density iy{z) is often approximated as a weighted combination of Dirac delta functions, 
as in 

cr(z) = ^ cr„(5(z„ - z), (6) 
n 

in which case 

VFr'(z)«^a„log(z„-z). (7) 

n 

Using a double-layer distribution instead of the single-layer distribution in ^ does not 
yield a second-kind integral equation, as it does in the classical formulations, but this is 
nevertheless a viable option 3) l^- The approach in this work is similar, in that we also 



use dipole potentials as basis functions. However, the source points are not restricted 
to lie on any given curve, but may be placed anywhere outside $7. Assuming that W{z) 
is also holomorphic for z G F, it can always be expressed as a series of complex dipole 
potentials, 

00 

W{z) = Y^^^, Vzen, (8) 

2„ Z 
n=l " 

where the z'^ are referred to as the poles of W{z). 

5 



According to a theorem by Wolff [40| (see also [4l|, Ch.l, Theorem 13]), given the 
appropriate choice of poles, z'^ ^ Q, ^ converges uniformly on any closed set in fl. In 
this sense, we may say that the set of single-pole functions is a complete set, which can 
be used to approximate any function W{z) which is holomorphic for z € fl. This theorem 
does not cover the case of W{z) having a branch-point on F. Nevertheless, in SectionjTU 
we show a solution with error < 10"^" even for such a case. The convergence theory of 
rational function approximations to functions with branch-points is rather involved, and 



we therefore refer the interested reader to a couple of relevant works [42|, |43|. 

Our algorithm is based on approximating the complex potential using N poles plus 
a constant, or equivalently as an N^^ order rational function 



a„ Pn{z) 



^^z^-z Qn{z) 

where Pn and Qn are complex-coefficient polynomials of degree N. To determine 
oq, oi, . . . , a^r, or equivalently, the coefficients of Pm{z)^ we start with an initial guess for 
the poles, and associate a complex basis function with each pole, 

= -7^, (10) 

z„ z 

with real and imaginary parts, 

0«(z) = Re[0^(^)], ^liz)^Im[^^iz)]. (11) 

Note that (j)n{z) and — ^^(z) are real-valued basis functions that correspond to x- and 
y-directed dipole potentials. In the special case of the constant term, the complex basis 
function becomes a real- valued basis function, 



(12) 



Defining U{z) = Re[M^(z)] and V{z) = lm[W{z)], and using ([T(I|-p^. 

N N 



U{z) « Uiz) = ^ a^cb^iz) - J2 <4'iiz), (13) 



n—O n— 1 

and its harmonic conjugate, 

N N 



n=l n—1 

where the basis function coefficients = Re(a„) and a!^ = Im(a„) are determined by 
minimizing an error metric on the boundary T. More precisely, 

a=[a^,af + ia\,...,a%+ia\!]'^ (15) 

is determined by solving the minimization problem 

U^f =mm(u-f,U-f)\ (16) 
r a \ / r 
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where the norm, || • ||r, is induced by a Hermitian inner-product, (•, •)p, the definition of 
which is deferred to Section [3] The minimization in (jl6l) is equivalent to requiring that 
the residual, f — U,he orthogonal to the span of the basis functions, or equivalently, 

(0^,/-^)^ = O, n = 0,l,...,7V, (17a) 
{^iJ-u)^^0, n = l,...,N. (17b) 

Expanding the inner products and using leads to the normal equations. 





$Ri' 










$IR 


$11 




ai 







(18) 



where ^± = (O^r and (fp„ = (/, 0^)p. 

Finding the basis function coefficients by least-squares minimization is common in 
desingularized methods. What is less common is our use of dipole potentials as basis 
functions, and also our use of a continuous error metric, as in (jl6p . In Section |3] we 
will define a Hermitian inner product, and show that under mild conditions, enforcing 
the residual orthogonality (1171) is equivalent to enforcing the boundary conditions at 
pole-location-related quadrature points. 



2.2. Pole relocation 



In the next step of the algorithm, the poles are moved to improve the approximation 
of W{z). This is done by fitting a rational function to the complex function 



W{z) = f{z) + iV{z), zeT 



(19) 



where V{z) in (fT4)) is based on the current set of poles. To obtain a "better" set of poles, 
we seek a rational approximation to W{z) in the form P{z)/Q{z) where both P{z) and 
Q{z) are polynomials of degree < N, and then use the zeros of Q{z) as the "better" 
poles. 

This rational fitting is done by minimizing a norm induced by the same inner-product 
used in (TTf)) . Ideally, we would want to minimize ||P/(5 — Vl^||r, but this is a nonlinear 
and non-convex problem. So, instead, we follow the commonly-used strategy of first 
scaling the residual P/Q — W by Q, and then correcting this scaling by dividing by 



N 



(20) 



n=l 



which is an approximation for Q{z) based on the current set of poles. The minimization 
problem is then 

P-QW 



mm 



(21) 



where P and Q are represented with polynomial bases, denoted p„ and Qn, n = 0,1, . . . , N , 
respectively. To avoid the trivial solution, Q{z) is assumed to be monic in the g„ basis. 
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The minimization in (j2ip leads to the orthogonahty conditions, 



\Q Q \ Q Q /r 

The choice of polynomial bases determines the conditioning of ([22]) , and, as the monomial 



basis is a particularly poor choice, other choices have been proposed [44|-|46 1 



The above approach to rational fitting has a long history, dating back to Levy [47 



and Sanathanan and Koerner |48j . and its application to real functions is much older. 
Usually, however, the minimization problem in (j2ip corresponds to fitting a rational 
function at a finite set of points on the imaginary axis or the unit circle. For that more 
common case, many implementations have been suggested, including vector fitting (VF 
iterated rational fitting (IRF) 0, \^ , and recently, robust rational least-squares 
We use VF as the default method, but in some cases it becomes poorly conditioned and 
we switch to a variant of IRF, the details of which are given in [Appendix A| 

The main difference between VF and IRF is the choice of basis polynomials. In VF, 
the choice is Pn{z) = (?n(^) = ^niz), where 

{Q{z) if n = 

For n = 1,2, . . . , N , the €„(z) are, up to a scaling factor, the Lagrange polynomials for 
interpolating at the z'^. As these polynomials are zero at all poles but one, they usually 
lead to a fairly well-conditioned linear system. However, if two or more poles are very 
close, conditioning may deteriorate. If this happens, we switch to using IRF, in which 
two polynomial bases are derived. For the Pn{z), we use polynomials orthonormal with 
respect to the inner-product 

{u,v)^p^(^,^\ , (24) 



Q Q 



r 



and for the qn{z) we use polynomials orthonormal with respect to the inner product 

, , /Wu Wv\ 

<"'">™ = \T't)/ '''' 

As a result of this choice, the linear system to be solved is made-up of two blocks of 
orthonormal columns, one for P and one for Q (see [Appendix A] ) . 

Once Q{z) is obtained, its roots are computed and they become the poles for the 
next iteration. Using these poles, U{z), the approximate solution of ([2]), is recomputed 
by solving (fTS|) and then using ([T^ . Similarly, using (fn| . V{z) is recomputed, W{z) 
in (fT9|) is updated, and the poles are relocated again. 

When computing the roots of Q{z) we may encounter roots inside fl. Although we 
do not discard these poles, their corresponding basis functions are not used to form U{z) 
and V{z). Hence, they are tacitly excluded from (IT3l) - (IT8l) . As the algorithm progresses, 
these poles usually migrate outside of Vl. Any poles that remain in when the algorithm 
stops are discarded. 



2.3. Initial guess for the poles 

The proposed algorithm is found to be quite insensitive to the initial pole locations. 
We set them using the following heuristic, but even random initial locations seem to do 
fine. The initial guess for the is obtained by solving ((22)) assuming W{z) — f{z), 
which is purely real, Q{z) = 1, and then allowing P and Q to be of degree 2N . Out of 
the 2N zeros of Q we pick TV zeros that are outside f2, and, if there are not enough of 
these, we supplement them with zeros inside Q, so that we have N zeros in total. These 
N zeros are used as the poles for the first iteration. The justification for using 2N poles 
along with a purely real W{z) is based on considering the fitting problem when i7 is a 
disk. For a disk, a purely real W{z) can be constructed with N poles inside VI and N 
poles outside il. Using only the poles outside, we obtain a function that is holomorphic 
in the disk and its real part fits f{z) on T. 

2.4. Algorithm summary 

The main steps of the algorithm are as follows: 



Algorithm 1 Solve simply-connected Dirichlet problem 

1: Set the 2^ to the initial guess, as explained in Section [2?3l 

2: while stopping criterion not met do 

3: Solve ([IH]) for the and a\. Use only ^ Q.. 

4: Use the and ajj to form V according to (fT4|) and W according to (|T9)) . 
5: Form Q according to (|20)l . use ah z'^. 
6: Solve dig) for P and Q. 
7: = roots(Q). 

8: end while 



3. Inner product definition and quadrature rules 



In each iteration of Algorithm 1, two norm minimization problems must be solved; 
one to determine an approximate complex potential given the poles, (|16p . and one to 
update the poles using the approximate complex potential, (pij) . Both problems can 
be formulated as follows: given a set basis functions tpi, ip2, .. - ipN^ defined on F, and 
a function to approximate, g, also defined on F, we wish to determine the vector of 
coefficients c that solves the norm minimization problem 



(26) 



The unique solution to this problem is given by 



(27) 



where = (V'm,'0n)r ^^'^ iS'4>)m = (5, V'm)r- It can be easily shown that 'S' is Hermi- 
tian positive definite and hence invertible, provided that the ipn are linearly independent 
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(that this is so for the 0^ and cfi^ follows from a trivial modification of Theorem 4 in 
where it was proved in a 3D context). The inner product used in this work is 



u{z{s))v* {zis))X{s)ds, (28) 



where A(s) denotes a positive weight function, and z(s) is a parametrization of T. In 
view of this definition, the inner-product integrals in ^ and may be evaluated using 
a quadrature rule. If (pSj) can be evaluated with a single quadrature rule for all u's and 
w's occurring in inner-product integrals in ^ and g^, it is simpler and more numerically 
stable not to use (P7)l explicitly, but instead to sample y/ X{s)g{s) and the y/ X{s)ipnis) 
at the quadrature nodes, and to set-up an overdetermined linear system for c. 





Cl 











(29) 



where si, S2, ■ • ■ , sk^ are the quadrature nodes and Ai, A2, . . . , are the quadrature 
weights. Equation (p9)) may also be written as 

A^$c«A5g, (30) 

where 

A = diag(Ai,A2,...,A/f^) 



g = [5(si),5(s2), • • -.gisK^)] 

^ mn — '0n(5m)- 

The least-squares solution to ([50)1 is given analytically by 



(31) 



(32) 



where superscript H denotes Hermitian transpose. It is more numerically stable, however, 
to obtain this solution by QR decomposition with column pivoting applied to (|29p . If the 
quadrature rule is exact, then the solution given in ([5^ also solves (pS)) (see, e.g., [22, 
Ch. 4]). 



3.1. Quadrature rules for (|16p 

The basis functions used to solve (ITB)) . restricted to F and seen as a function of the 
parameter s, are the real and imaginary parts of 



1 



(33) 



Clearly, constructing quadrature rules that will integrate inner products of these functions 
for arbitrary z{s) is not practical. However, the parametrization of z{s) is ours to choose, 
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and we choose it to be a piecewise rational function, i.e., 



zis) 



'ri(s) if < s < hi 
r2{s) \i hi < s < h2 

jAiis) if hM-i < s < 1 



(34) 



where the r^is) are rational functions in s. This choice of parametrization is quite 
general, as piecewise rational functions can be used to represent any geometry of prac- 
tical interest, and they are easily generated from T using readily available libraries (we 
use both Matlab's spline toolbox and the Chebfun system 5l|). With this choice of 
parametrization, the (/>^(s) are themselves piecewise-rational functions, each piece hav- 
ing known poles in the complex s plane. Quadrature rules for products of rational func- 
tions with known poles are readily available, and this is the motivation for our choice of 
parametrizatiorQ. However, the actual basis functions used in (1161) are (/'JJ (s) and (/'Ji(s), 
and these are not piecewise-rational functions of s, as they depend on both s and s* . To 
overcome this difficulty, we use the Schwartz reflection principle. We write the real basis 
functions as 



1 


1 




Wn^riis)]* 


1 


1 


<- 7-2(5) 


[< -^2(5)]* 


1 


1 


4-^m(s) 


K - rM{s)]* 



ii ho < s < hi 
if /ii < s < /l2 

if hM-i < s <hM 



(35) 



where ho = 0, hM = 1- Note that the expression for cj)}^ is similar. Although (j>^ and 
(^J^ are functions of both s and s* , for real s they can be written as piecewise rational 
functions in s alone. These rational functions can be obtained by replacing [z'^ — rm(s)]* 
by [z'^ — r„i(s*)]* everywhere in psp. and similarly for (pl^. We then have. 



n 

1=1 



for s e [hm,hm+l), 



m 



0,1,...,M-1, (36) 



(s - aimn) (S ■ 



(=1 



where L„i is the degree of the numerator or denominator of r„i{s), whichever is higher. 
In ([55)1 . airan is the l^^ pole associated with the n"' basis function, restricted to the m}^ 
section of the boundary, i.e., it is the l^'^^ zero of z!^ — rm(s). 



^Another advantage of this parametrization is that determining whether a pole is inside or outside 
of f2 is easy (see [Appendix B] ). 
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Quadrature rules that are exact for the inner products of the basis functions are 
readily obtained using rational Gauss-Chebyshev quadrature rules as building blocks. 
These rules integrate exactly integrals of the form 



R{s)T*{s*)^==ds, (37) 
V 1 — s 

where R{s) and T[s) are rational functions that have poles at known locations (and need 
not be proper) [52|. Once the aimn are determined, a quadrature rule for each section of 
the boundary is calculated, in 0{NLm) time, using the Matlab toolbox Rcheb [HI]- The 
rules are then concatenated to obtain a single rule for the entire boundary. 

Using these quadrature rules implies a specific choice for the inner-product weight 
function, denoted A(s) in ([25)) . namely. 



A(,s) ^ < 



Ai(s) if ho < s < hi 

A2(s) if /ii < s < ft,2 

. , (38) 

Am(s) if h-M-i < s <hM 



where 



Xm[S) = ^ ■ (39) 

{hm+1 — ^m)^ — (2s — hm+1 — /ijn) 



This choice of A(s) is dictated primarily by the ease with which rational Gauss-Chebyshev 
quadrature rules may be computed. However, it also has the advantage that when 
the boundary consists of a single rational-function section, using the Chebyshev weight 
function implies that the approximation will be close to optimal in the max-norm sense, 
provided the boundary data is smooth enough (Hil . Issj . When the boundary is composed 
of more than one section, this is no longer so, but the Chebyshev weight function still 
helps suppress Runge's phenomenon near s = and s = 1. Lastly, this choice of weight 
function makes our error metric dependent on the parametrization z(s). To avoid this, 
one could choose A(s) = \dz/ds\, which would make the inner product a line integral and 
hence invariant to the choice of parametrization, but this choice would complicate the 
construction of the quadrature rules. A compromise is to make the weight of each section 
of the boundary proportional to its length, and this can be done with the chosen weight 
function by making the size of the interval [hm, hm + 1] proportional to the length of the 
m}^ section of the boundary. 

We still need to ensure that the inner-products between basis functions and the 
boundary data / {z{s)) are calculated accurately. Assuming the boundary data is smooth 
on each section of the boundary, this can be accomplished by making the quadrature 
rules exact for improper rational functions with a sufficiently high numerator polynomial 
degree. In Rcheb, this is done by adding poles at infinity when generating the quadrature 
rules. 



3.2. Quadrature rules for problem (1211) 

The inner products for this second norm minimization problem are given in (j22p . All 
the integrands in this case can be written as rational functions of s, with poles due to 
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the zeros of Q {z{s)) as well as the poles of W {z{s)). Recalling the definition of W {z{s)) 
given in (jl9p . its real part is just the boundary data, and hence is assumed piecewise 
smooth, while the imaginary part has the same poles as those used for p^ . To these 
poles we must add all zeros of Q (z(s)), which are the zeros of all the z'^ — rm[s), including 
any poles z'^ G fi. Although the quadrature rules are calculated quickly by Rcheb, it is 
usually more economical to call it only once, with the poles for (PT|) and then to solve 
both equations with the same set of quadrature points. 

The proposed scheme for evaluating inner-product integrals with quadrature yields 
linear matrix equations very much like those of conventional desingularized methods, 
except that quadrature dictates the location and weight of each testing point. This 
scheme takes into account the location of the sources, the shape of the boundary, and 
the variation of the boundary data, and minimizes a continuous error norm. 

The total number of testing points used, K , is given by 

M 

K = M{N^ + 1) + (3iVout + M„) (40) 

m— 1 

with the notation 

M number of rational-function sections of the boundary 

Naa number of infinite poles added 

iVout number of poles outside of 

TVin number of poles inside of 57 

Lm degree of numerator or denominator of r„i(s), whichever higher. 

As the scheme is based on Gaussian quadrature, this number is in a sense optimal for 
computing with the continuous error norm, and the computation of the testing points is 
done with a fast and accurate library routine. We have found that, in all but the sim- 
plest cases, many uniformly distributed testing points are needed to match the accuracy 
obtained by the Gaussian quadrature scheme, so the extra computational burden is gen- 
erally worthwhile. However, the representation of the boundary with piecewise rational 
functions should be made as simple as possible, i.e., with the fewest sections and lowest 
orders. If the boundary is broken into many sections, then taking into account all the 
poles when computing the quadrature rule for each section can be excessive. It is better 
to take into account only poles that are close to a given section. The rest of the poles 
can be dealt with more efficiently by a small number of poles at infinity. We refer to this 
process as pole pruning. 

3.2.1. Pole pruning 

Adding N^a poles at infinity makes the quadrature rule exact for integrands that 
are products of a rational function and a polynomial of degree 2Noo — 1. Hence, we 
can omit poles that are distant enough that their corresponding potentials are smooth. 
To determine which poles are far enough, we need to estimate the error committed in 
approximating a rational function by a polynomial on the integration interval. As a 



13 



prototypical rational function we take a third-order pole function, 



1 



(41) 



as this is the highest pole order that will appear in the integrands, and estimate the error 
committed in approximating it by a polynomial, as a function of the pole sq. The max- 
norm error of the polynomial that interpolates at Chebyshev nodes is almost minimal 
and it is the same for all sq on the ellipse 



s(^) 



pe 



G [0,2^], 



(42) 



where p > 1 is the so-called elliptical radius (see, e.g., [54)|). For any Sp ^ [— Ij 1] we can 
calculate the elliptical radius corresponding to the ellipse which passes through it. 



p(so) 



So + 



1 



So 



1 



(43) 



Although there are bounds for the error in approximating a function such as i?so (s) by a 
polynomial, we found it more practical to tabulate the degrees of polynomials that yield 
approximation errors less than e, for relevant values of e and p{sq) (Fig. [3]). To obtain 



400 




logio(p-l) 



Figure 3: Polynomial degree that yields approximations of Rsg{s) with errors less than 
e, as a function of p. For p « 10, errors under lO^^'' can be obtained with polynomial 
degrees under 20. 



these polynomial degrees, we again used the Chebfun system, which truncates a Cheby- 
shev expansion once the expansion coefficients drop sufficiently below e. The tabulation 
is then interpolated or extrapolated as necessary so that we can rapidly calculate the 
optimal p beyond which poles can be discarded and replaced by poles at infinity, without 
incurring errors much larger than e. As shown in Fig. |31 for all but the most distant 
poles, the choice of e influences the required polynomial degree appreciably. Hence, if 
only moderate accuracy is sought, significant reduction may be gained by using a larger 
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e. To take advantage of this, we set e according to 

e = min(10"'*,0.01A£;), (44) 

where IS.E is the normahzed error obtained when solving (jl7p in the previous iteration (a 
more precise definition is given below in (1551) ). Thus, if the solution is not very accurate, 
fewer quadrature nodes are used, but their number is increased if the algorithm is allowed 
to iterate and higher accuracy is achieved. In the very first iterations, /S.E can be large, 
so we bound epsilon from above by 10"''. Of course, this rule is heuristic and has been 
derived from numerical experiments in which it appeared to do well, but it is by no means 
optimal and a different rule could perform better. 



3.2.2. Froissart doublets 

Admittedly, the quadrature scheme just described is somewhat complicated, but we 
have found this complication to be worthwhile. The reason is related to the occurrence of 
pole- zero pairs, called Froissart doublets [55-57|, which almost coincide with the bound- 
ary. These pole-zero pairs are a numerical artifact, and they do not contribute to the 
fidelity of the fit, but instead they fit errors in the data. In our case, the errors in the 
imaginary part of the data are quite large in the first iterations, so these artifacts occur 
often. If uniformly distributed testing points are used, the doublets frequently occur 
closer to the boundary than the separation between testing points. When the poles of 
these doublets are used to solve (ITBl) , the resulting potential is inaccurate between testing 
points. In turn, the estimate of Im[VF(z)] is inaccurate, and the improvement in accuracy 
from one iteration to the next may stall. One could try to filter these poles in various 
ways, for example, based on their small residues and distances from the boundary, but 
this would limit the obtainable accuracy. By using the proposed quadrature scheme, all 
the poles outside of £7, no matter how close they are to its boundary, can be safely used 
to solve P^ . 

In Fig. Hh., an artifact caused by a Froissart doublet can be observed near the top 
of r. The pole of the doublet is extremely close to the boundary (the distance is about 
10^^ of the perimeter of F), which explains the appearance of the artifact despite having 
used 2000 uniformly distributed testing points. When the quadrature scheme is used, 
with just 203 points, the artifact disappears (Fig. HJj). 



4. Dirichlet problem, multiply-connected domain 

Other Laplace BVPs can also be solved with the proposed scheme, after appropriate 
adaptations. As an example of a Dirichlet problem in a multiply-connected domain 
we consider the exterior Dirichlet problem; multiply-connected interior problems can 
be treated similarly. In the exterior Dirichlet problem, the potential U{x,y) obeys the 
Laplace equation everywhere outside of a set of disjoint domains f2i, r22, • ■ • , bounded 
by contours Fi, . . . , Fj, respectively. As before, Dirichlet boundary conditions are 
imposed, 

C^lonr, i = l,2,...J (45) 

and in addition, the potential must be bounded at infinity. Of course, the poles must 
now be restricted to lie inside the VLj. Also, the complex potential W{z) is no longer 



15 



(a) 



(b) 



Figure 4: Equipotential lines for a monopole (marked by a red dot) inside an L-shaped 
region. In (a), 2000 points were uniformly distributed on F, in (b), 203 points were 
distributed by the quadrature scheme. 



single- valued, but can he written as the sum of J logarithmic terms and a single- valued 
holomorphic function W{z) (see, e.g., ^3, §29]). We have, 

J 



W{z) = W{z) +J2 A, log (z-zf), 

3 = 1 



zfen,, 



(46) 



where the Aj may be assumed real, and the zj^ can be anywhere inside flj, preferably 
as far from its boundary as possible. If the Aj were known, the logarithmic terms 
could be subtracted from the boundary jiata, and the algorithm for the simply-connected 
Dirichlet problem could be applied to W{z). Our approach is therefore to estimate the 
Aj iteratively, along with the imaginary part of W{z). On each iteration, we solve 



Tl=l j = l 



z — z. 



71=0 



(47) 



with /(z) = jj(z) for z € Fj, 
J 

subject to E =0, 

for the aJJ, the a^, and the Aj. The condition that the sum of Aj's is zero, needed to 
ensure the solution is bounded at infinity, can be enforced by eliminating one of the Aj , 
say, Aj = — X^fji ^i- Instead of (H71) we then have. 



Naut Afout J-1 

E «n'/>^(^) - E ""n^ni^) + E 
n— n—1 i — 1 



(48) 



As before, we solve (1481) by sampling at quadrature nodes, weighting the equation by 
the square roots of the quadrature weights, and solving the linear system in the least- 
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squares sense. In contrast to the previous case, however, the inner products involving the 
logarithmic basis functions will not be integrated exactly. We must therefore add poles at 
infinity when deriving the quadrature rules for each section of the boundary. Fortunately, 
the logarithmic singularity is weak, so when the are not too close to the boundary, 
the restriction of the logarithmic term to the boundary curve can be approximated well 
by piecewise polynomials. The poles at infinity added by the pole-pruning process of 
Section [3. 2. II usuallv sufHce, although extra poles at infinity can be added if necessary. 

Once (|48p has been solved, a rational function is fitted to an estimate of the boundary 
values of Wi^z), 



J-i 

Ty(z)«/(z)- V^.log 



(49) 



The real part is estimated as ](z)^ the known values of Re[W^(z)] on the boundaries, with 
the real part of the logarithmic terms subtracted. The estimate for Im[Vl^(z)] is V"(z), 
given by (O, and it remains as it was for the simply-connected problem when W(z) and 
W{z) were the same. 



5. Exterior Neumann problem 

Now consider the case when the normal derivative of the potential on the boundary 
is prescribed. For brevity, we discuss only the exterior problem, which can be formulated 
as 



V2c/(a;,y) = 0, for all (x, y) ^ (J 



dn 



(50) 



on 



lim y)\z\ < oo, 

|z|— >CXD 



^y, 



where the flj are again simple disjoint domains bounded by the Tj. This problem is 
uniquely solvable if and only if 



(51) 



The approach we adopt for this problem follows [33, §32]. As in the exterior Dirichlet 
problem, we seek the solution as the real part of a function W{z) that is holoniorphic 
outside the flj, and which can be written as a sum of a single- valued part, W{z), and J 
logarithmic terms. In this case, however, the logarithmic terms are known in advance: 



(52) 



17 



and they can be subtracted from the boundary data. After this subtraction, 



^ J 

f J (z) ^ Mz) A, log \z-z^\, (53) 

can be used to write a Neumann boundary condition for Re W{z) . Denoting W{z) = 
U{z) 4- iV{z), this condition reads, 



dU 
dn 



h- (54) 



If W{z) is approximated by a weighted combination of dipole potentials, it becomes 
inconvenient to enforce (IMl) directly, as the differentiation increases the order of the 
poles. Instead, we integrate ((54)) with respect to arc-length and use the Cauchy-Riemann 
equations to obtain. 



V{z) = / f,{z)\dz\ + V, EE F,{z) + (55) 

Zj 

where the integration is along Tj, from an arbitrary point Zj to the point z, and the 
Vj are constants yet to be determined. Finding a harmonic function V that satisfies 



the boundary condition (|55p is known as a modified Dirichlet [37[ or Muskhelishvili f5_8| 
problem, and given the condition at infinity in (j50p . it is solvable only for a unique set 
of Vj constants. These can be determined along with an approximate solution for V , by 
solving 

^aJj0;,(z) + a5>^(z)-i;,«F, (z), for z € F„ (56) 

n=l 

for the a^, the a^, and the Vj. As before, (|55)) is solved in the least-squares sense after 
sampling the basis functions at quadrature nodes and weighting the equations by the 
square roots of quadrature weights. Once ([56)) is solved, the rational fitting algorithm is 
applied to an estimate of the boundary values of W{z) given by. 



_ A^out 

iy(z)« ^aJ,(/)J,(z) + a^0^(z) + z[F,(z) + t;,], for z G F,. (57) 

n=l 

As expected, in the Neumann problem, it is the imaginary part of the complex potential 
that is prescribed by the boundary data (up to J additive constants), while its real part 
and the J constants are estimated iteratively. Note that the ti = term corresponding 
to the constant basis function is omitted from the sum in ()57p so as to conform with the 
boundary condition at infinity in ()50p . 

The formulation just described differs from one that is commonly used for Neumann 
problems in integral-equation formulations (e.g., in [38|). These formulations typically 
use a single-layer potential so that enforcing the boundary condition of ([50)) leads to a 
second-kind Fredholm equation. We could have formulated our solution similarly, i.e., by 
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representing the potential as a weighted combination of monopole potentials. Enforcing 
the boundary condition would then lead to a weighted combination of dipole potentials 
approximating the normal derivative of U. The chosen formulation does have some 
advantages, however. First, the functions to be approximated, the Fj{s), are smoother 
than the fj{s), so fewer basis functions may be needed for their approximation. Second, 
a common situation in Neumann problems is that the Aj of ([5^ are zero, so that the 
complex potential iy(z) is single valued. This occurs when the boundary is composed 
of a single contour, or when the Qj represent impenetrable objects. In these cases, it is 
convenient to approximate W{z) using single-valued basis functions, instead of having 
a superposition of non-single-valued functions that can be made single-valued by an 
appropriate choice of branch cuts, as would occur if monopoles were used. For instance, 
using the proposed formulation, it is easy to obtain equipotential contours and field lines 
as the level sets of Re [W^(z)] and Im without having to keep track of the proper 

Riemann-sheet on which to evaluate the logarithmic kernel. 

6. Large-scale problems 

The proposed algorithm for setting the poles and testing-points works well for moderately- 
sized geometries, but its computational cost is too high for larger ones. However, it is 
still possible to tackle large-scale problems by partitioning the geometry into smaller 
domains and running the algorithm for each of these domains separately, possibly in 
parallel. Once the poles and testing points are obtained for the entire geometry, we 
solve the full problem with the locally-optimized poles and testing points. If the cost 
of solving the full problem becomes prohibitive, it may be possible to use an iterative 
domain-decomposition approach, along the lines of [H^ , [g^] , or [6l| , which could also be 
used to improve the locations of the poles and testing points from one iteration to the 
next. 

We propose a scheme for geometries made up of a large number of simple closed 
curves. The scheme consists of partitioning the boundary curves into clusters so that 
each curve belongs to a single cluster, as shown in Fig.[S] Various algorithms for obtaining 
such a partition exist, and we use the simplest one, the fc-means clustering, or Lloyd's, 
algorithm [62!]. To keep the partitioning simple, each curve is entirely within one cluster 
and represented by its centroid, on which the clustering algorithm operates. Then, to 
each cluster we add a "buffer region," as shown in Fig.O so that curves close to the edge 
of the cluster have all their nearby interactions taken into account. The buffer region is 
obtained by adding to each curve all the curves with centroids closer than some radius, 
and this radius is chosen as twice the distance to the nearest centroid among all other 
curves. Poles that, after the local optimization, are found inside the buffer region are 
not used when solving the full problem. 

7. Numerical results 

The proposed method was implemented in Matlab, and results for a number of sample 
problems were obtained. This implementation, which is freely available for download [63| . 
was run on a 2GHz Intel Xeon machine with 8GB of RAM. The error metric adopted in 
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Figure 5: Boundary curves partitioned into three clusters with buffer regions (dehmited 
by dashed hues). 



this work is AE, 



AE = 



u 



(58) 



where f and u are vectors of samples at the quadrature points of the exact and ap- 
proximate boundary values, respectively. The || • ||„, norm includes multiplication by the 
quadrature weights, so that it is an approximation of the continuous norm. The residual, 
U — U, is itself a harmonic function throughout and hence is largest at the boundary. 
We also used the maximum-norm error, 



max 

zer 



/(^) - U{z) 



ll/WllL^(r)/l|l| 



lL2(r)/l|i||L2(r) 

where ||l||i2(-p-) is just the length of F. In practice, we estimated AE„ 
and ?7 on a large number of points on the boundary. 



(59) 

by sampling / 



7.1. Interior Dirichlet Problem 

Problems of increasing difficulty were considered. Where applicable, results were 
compared to results obtained with the standard Nystrom method, implemented along 
the lines of [sSj]. However, we have not compared with more sophisticated Nystrom 
schemes, such as those of [l|, [Ij H 0, which have been recently proposed for handling 
more "difficult" geometries and excitations. 

7.1.1. A few well-separated poles 

We begin with the simplest problems, in which the exact solution is known in advance 
to have just a few, well-separated, poles, i.e.. 



f{z) - Re 



-No 

E 



-1 "fe 



- z 
20 



(60) 



Although this excitation suits our basis functions perfectly, the solution obtained is not 
always very accurate even if iV = A^o- If the estimate of Im [14^(2)] is inaccurate enough, 
the poles may be relocated to poor positions, and using these poles will not lead to a 
more accurate estimate of Im [iy(z)]. This problem can be alleviated by using N > Nq, 
though the number of poles that must be added beyond A'o depends on the location of 
the poles of the solution and also on the shape of the boundary. This dependence is 
explored below. 

A plot of error vs. iteration count is shown in Fig. [51 for an ellipse with unit semi- 
major axis and an aspect ratio of two, and boundary data as in (|60|) with A'o = 6, = 
k/No, and z'f, — 1.5 exp [i7r(2fc + l)/6]. When N = Nq, the error does not decrease 
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Figure 6; Error vs. iteration for an ellipse and 6-pole boundary data. One of the poles 
remains inside the ellipse, both for A^ = 6 and N — 7. 



much below 0.1%, but when A^ = Aq + 1, the error drops exponentially, approaching the 
machine cpsilon after 45 iterations. The reason for this dramatic difference in accuracy is 
that when A — Nq, a single pole remains inside ft and the estimate of Im[iy(z)] obtained 
with the remaining five poles is not accurate enough to make it migrate outside of Q. 
Consequently, the improvement in accuracy stalls. In contrast, when A^ — Nq + 1, the 
extra pole allows for a better approximation of the imaginary part and six poles appear 
outside of Q. They are then relocated until the error approaches the machine epsilon. 

A natural possibility is to increase A^ adaptively, whenever the improvement in ac- 
curacy stalls. For example, one can calculate (Tk{^E), the standard deviation of AE in 
the previous K iterations, and add a pole whenever 



CT/f (AS) < fiKiAE)e, 



(61) 



where fj.K{AE) is the AT-average of AE and e is a small number. A more sophisticated 
strategy might be found along the lines of [H^ , which could also be used to decrease the 
number of poles if necessary. We have implemented the simple strategy with A' = 5 and 
e = 0.5%, which appears to yield a practical tradeoff between the conflicting goals of 
only adding poles when necessary, and keeping the number of iterations moderate. This 
adaptive scheme, which makes the algorithm fully automatic, is available as part of the 



downloadable package 



However, to keep the results section concise, we only show 



results for the nonadaptive version. 
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We consider next boundary curves of the form, 

z{s) = [1/7 + (1/7 - 1) cos(27ri/s)] e^'^''", s e [0, 1] (62) 

with 7 = 1.75, and = 2,4, 6. A plot of the ly = 6 curve is shown in Fig. [71 together 
with a set of ten points distributed on a circle of radius R, which represent the poles of 
the boundary data. This radius was set to 1.1, 1.5, and 3, while the number of poles was 
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Figure 7: Boundary curve and poles of boundary data used in Table [TJ 

kept to ten. The algorithm was run for a maximum of 100 iterations and we determined 
the smallest N required so that AE < lO^^"' when the algorithm stops. The results of 
this experiment are summarized in Table [T] where this value of N is given, together with 
AE'max, estimated by evaluating the approximate solution at 20000 boundary points. In 



Table 1: Results for boundary curves as in Fig. [Tj 
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all cases, Ai^max is of the same order of magnitude as the tolerance for AE, although 
in a few cases it is slightly larger than lO-i"*. We see that, in contrast to conventional 
desingularized- and integral-equation- methods, the proposed method works best when 
the poles of / are close to the boundary, adding only a few spurious poles during the 
solution. Evidently, when the poles are further away their locations are more difficult to 
ascertain. We note that v has only minor influence on N. 
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7.1.2. Almost- singular boundary data 

When the boundary data is singular, or almost so, the accuracy of most numerical 
methods suffers. A possible remedy is to add basis functions that can approximate 
the singular boundary data to the set of basis functions being used (for an example in 
the context of desingularized methods, see [S^). In the proposed method this happens 
automatically, as the basis functions are chosen from a large over-complete set. In the 
next example, the boundary is an ellipse with unit semi-minor axis, and an aspect ratio 
of two, and the boundary data is. 



f(z) = Re <^ exp 



1 



1.01 



z er. 



(63) 



This implies that the exact complex potential has an essential singularity very close to 
the boundary. Nevertheless, errors on the order of Ai?,nax — 10^^^ are obtained with no 
more than 25 poles (Fig. [5]). These errors could not be reduced further, and in fact, as 
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Figure 8: Almost-singular boundary data. Comparison of the error vs. number of basis 
functions {2N) obtained with the proposed method and the Nystrom method. 

more poles are added the error increases slightly. Compared to the proposed algorithm, 
the Nystrom method yields very large errors, and its convergence rate is very slow. For 
200 collocation points the Nystrom method yields Aii^max — 0.38, and, outside the range 
of the figure, 1000 collocation points still yield only /S.E — 0.05. 

Of course, for the same number of basis functions, the computation time of the 
Nystrom method is shorter by a factor related to the number of pole relocation iterations. 
Even so, the Nystrom method is more time-consuming in this example for any reasonable 
accuracy. For example, for /S.E — 0.05, the Nystrom method takes about 6s, while the 
proposed method takes 2.5s and it uses five poles. For higher accuracies the differences 
increase and it quickly becomes impractical to use a direct solver with the Nystrom 
method. If the boundary data is smoother, for example, if the essential singularity is 
moved to z = 1.1, the Nystrom method can be faster. It requires 0.9s and N = 400 for 
/S.E = 10^^, whereas the proposed method requires 2.5s and five poles to achieve the 
same accuracy. 

At each point on the error curves in Fig.|H]we calculated the error of an A^-pole solution 
and a 2A^-point Nystrom solution, so as to match the number of basis functions. The 
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maximum error for both the Nystrom method and the proposed method was calculated by 
evaluating the approximate potential on 8iV boundary points. To obtain this potential 
in the Nystrom method we used a 8A^ x 8A^ collocation matrix, which multiplies the 
double-layer density interpolated trigonometrically for its values at the 8A^ points. 

7.1.3. Almost-sharp corners 

The influence of the shape of the boundary in the previous examples is limited to 
determining where the solution is sampled; it does not influence the location of poles of 
the solution as these are set by the boundary data. In the next example, we consider 
a Green's-function type problem, where a boundary curve of form (j62p with = 2 is 
excited by a monopole potential, 

/(s) =log|z(s)|. (64) 

As the parameter 7 — > 2 in (j62p the problem becomes more difficult, because the minimal 
radius of curvature decreases, and also the sharp tips of the boundary approach the 
monopole and each other. A solution of this problem with 7 = 1.9, obtained with 35 
poles is shown in Fig. 1^1 The inhomogeneous distribution of poles and testing points is 




Figure 9: Equipotential contours (logarithmically spaced) and field lines for a monopole 
inside the curve: z{s) = [1/1.9 + (1/1.9 — 1) cos(47rs)] e^'^^'^ . The poles of the solution are 
marked by green circles, with centers colored according to the magnitude of the residues 
(darker means larger residues). The testing points are marked by blue dots. 

clearly necessary in order to approximate the rapidly varying field where the radius of 
curvature is small. To calculate the testing points, the boundary was parameterized by 
a single polynomial in s, of degree 34, which approximates /(s) to machine precision. 
Equipotential contours and field lines are also shown in Fig. [HI These were obtained as 
the level sets of Re [VF(z)] and Im [14^(2:)] by use of Matlab's 'contour' command. The 
evaluation of W{z) throughout is trivial, as it only requires the summation of the 
partial-fraction expansion for W{z). This may be contrasted with the more elaborate 
schemes needed to evaluate integral equation solutions throughout the domain, especially 
close to the boundary (see, e.g., [H). 

Plots of Ai?inax VS. N are shown in Fig.[TU]for curves with 7 = 1.5, 1.75, 1.9, together 
with results obtained with the Nystrom method. The points for the Nystrom method 
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were distributed uniformly in the parameter s of (I62p . implying a higher density near 
the origin, which helps in resolving the rapid variation in the double-layer density there. 
The proposed algorithm was allowed to run for a maximum of 100 iterations, though the 
poles stopped moving after roughly 50 iterations. The rate of convergence with respect 
to the number of poles is approximately exponential, and the minimal error obtained 
is similar to that of the Nystrom method. It is also evident that the proposed method 
generally yields a more accurate solution per basis function than the Nystrom method, 
and that the difference between the two increases with 7. 

7.1.4. Perfectly sharp corners 

When the boundary has sharp corners, the potential may have singularities at the 
corners, making it impossible to represent with simple poles alone. It is still possible, 
however, to obtain very good approximate solutions with a moderate number of poles. 
Such an example is shown in Fig. [T|d, and another one is shown in Fig. 1111 As can be 
observed in both figures, the poles and testing points cluster near the reentrant corners, 
and appear to be approximating a branch-cut which exists in the analytical continuation 



of the solution (see [43| for a discussion of the curves on which poles cluster when ap- 
proximating algebraic functions by rational functions). The distribution of poles near the 
reentrant corner is examined in Fig. 1121 where the poles are plotted in descending order 
of distance from the corner, for a solution with 82 poles of the L-shaped boundary. The 
distances, which are normalized with respect to the distance of the exciting line-source 
from the corner, are seen to decrease at an almost constant exponential rate, with the 
closest poles clustered even more densely. The closest pole is placed at a normalized 
distance of about 10""'^^, indicating that poles must be located very precisely if high 
accuracy is to be achieved. Moreover, it seems inevitable that the testing points must 
also be adapted so that they are located with comparable precision. One might wonder 
whether the sharp turn in the equipotential lines near the corner is achieved with poles 
having large residues of differing phases. As shown in Fig. [T^l this is not the case. Al- 
though the residues have absolute values that span 19 orders of magnitude, the smallest 
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Figure 11: Equipotential contours (logarithmically spaced) and field lines for a monopole 
inside an L-shaped domain. 

ones correspond to the closest poles. 

We now examine the convergence of the algorithm with N (Fig. [T^ . As seen in 
Fig. [T^ . the dependence of the error on N is not as smoothly exponential as in the pre- 
vious examples. This more uneven error decrease is due to inaccuracies in the imaginary 
part that result in poles converging to locations inside il. For example, when = 20 
four of the poles remain inside fi, while when = 30 all poles move outside of $7. As a 
result, AE decreases by more than two orders of magnitude when going from A^ = 20 to 
N = 30. As shown in Fig. 113b . the condition number increases with A^ , but the norm of 
the vector of residues (normalized with respect to A^) decreases. This is similar to [25j, in 
which such behavior was observed when the sources enclosed all the singularities of the 
analytic continuation of the potential. It is remarkable that this occurs in the present 
example as the potential cannot be continued analytically beyond the corners. 

7.2. Exterior Dirichlet problem 

The exterior Dirichlet problem differs from the interior problem by the presence of 
logarithmic terms in the former. Also, the region in which the poles may be placed 
in the exterior problem is bounded, which probably makes their location easier. A 
typical difficulty in exterior problems occurs when the boundary is composed of two or 
more curves that nearly touch. An example of this sort consists of two circular curves 
with boundary values —1 on the lower curve and +1 on the upper one, and relative 
separations d/R = 0.01,0.1,0.5, as shown in Fig. [T31 This example has an analytic 
solution [gHI to which we may compare the results. Specifically, seen as a two-conductor 
electrostatic problem, a global parameter of interest is the capacitance, proportional to 
the magnitude of the logarithmic terms, and a local parameter of interest is the charge 
density, proportional to the normal derivative of the potential. As can be seen in Fig. 1141 
for both of these parameters, the proposed method requires far fewer degrees of freedom 
than the Nystrom method, and this difference increases as the conductors approach. 
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Figure 12: Distances to reentrant corner, normalized with respect to the distance of 
the exciting hne-source from the corner, and corresponding absolute value of residues, 
normalized with respect to the largest one. 



Also of note is that the proposed method is far less sensitive to the separation between 
conductors. To a good approximation, the N needed to achieve a given error scales like 
1/d when using the Nystrom method, but it scales like log(l/d) when using the proposed 
method. 

7.3. Large-scale problems: exterior Neumann example 

In the last example, we solve a Neumann problem of potential flow around impenetra- 
ble ellipses of random radius, eccentricity, and orientation, many of which are very close 
to one another. This geometry was generated by a simple modification of the method 
described in 0] for generating a distribution of random circles. The velocity field far 
from the ellipses is assumed uniform, of unit amplitude, and in the x direction. The 
corresponding BVP is then given by ((50)) with fj{z) — — cos[9j{z)] where Oj{z) is the 
angle between the outward-pointing normal to the j*^ curve and the x axis. 

We use the partitioning scheme of Section [B] and set the number of clusters so that 
each cluster has only a few ellipses. Thus, the run-time and memory requirements for 
each cluster are kept very moderate. For example, for a geometry with 100 ellipses, we 
used ten clusters, with the largest cluster having 13 ellipses (not including the ones added 
as buffers). For each cluster, we used three poles per ellipse (including the ones used 
as buffers), and stopped the algorithm if the error dropped below 1% or 100 iterations 
were exhausted. The algorithm was run in parallel on four cores of a desktop PC, using 
Matlab's parallel computing toolbox, and the run-time for determining the pole and 
testing point locations was 87s. We then assembled a linear system of size 4832 x 706 
with the poles and testing-points of all clusters and solved it by QR decomposition in 5s. 
The error of this solution was AE = 0.013, slightly higher than the error tolerance for 
each cluster. The largest geometry solved, which consisted of 300 ellipses, is shown in 
Fig. [13 together with computed flow lines. Here we used 35 clusters, and the run-time 
for determining the pole and testing points was 9 minutes. The resulting linear system 
was of size 14581 x 2150. Although fairly large, this system has, on average, 7.2 degrees 
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(a) 



(b) 



Figure 13: L-shapcd domain, monopole excitation. Despite singularities on the boundary, 
errors smaller than 10~^° can be obtained with no more than 100 poles (a). Also, 
although the condition number increases with N, the norm of the vector of residues 
(normalized with respect to N) decreases with N (b). 
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Figure 14: Two closely-spaced curves. Full lines, A^-pole solution. Dashed lines, 2N 
Nystrom solution. The accuracy of the proposed method is far less sensitive to the 
separation between circles than the Nystrom method. 



of freedom per curve, and it was solved in 2 minutes by QR decomposition. The error 
tolerance for each cluster was again AE < 0.01 and the error of the full solution was 
AE = 0.014. 



8. Summary 

In this paper, we describe an algorithm for solving 2D Laplace EVPs with a desin- 
gularized method. The algorithm is based on rational-function fitting techniques, and 
it uses rational Gauss-Chebyshev quadrature rules to locate testing points at which the 
boundary conditions are enforced in a least-squares sense. To use these rules, the bound- 
ary curves are parameterized by piecewise-rational functions, which can be used to ap- 
proximate quite arbitrary boundary curves. The algorithm can be viewed as an adap- 
tive desingularized method, in which the sources are not restricted to curves conformal 
with the boundary curves, but can be placed anywhere outside the EVP domain. This 
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Figure 15: Stream lines around 300 ellipses. Calculation using the partitioning scheme 
took 11 minutes on a single 2Ghz PC with a four-core processor. 



added flexibility allows the method to adapt to difficult boundary geometries and almost- 
singular boundary data, yielding solutions with errors approaching the machine epsilon 
for a relatively small number of degrees of freedom. Examples shown include boundaries 
with sharp and almost-sharp corners, boundary data that is derived from functions with 
essential singularities, and boundary curves that nearly touch. We also described a par- 
titioning scheme, in which pole optimization is limited to curves within a cluster, that 
can be used to tackle large-scale problems. 

As the proposed method relies on a correspondence between analytic functions in 
the complex plane and harmonic functions, and on properties of rational functions, it 
may be difficult to generalize to 3D Laplace problems. On the other hand, generalizing 
the method for 2D Helmholtz or Biharmonic BVPs appears more straightforward. The 
general solution of the Biharmonic equation can be expressed in terms of two functions 
that are analytic throughout the BVP domain, and rational functions could be used for 
their approximation. For Helmholz BVPs, it may be possible to exploit a connection be- 
tween the singularities of the analytic continuations of solution to Helmholtz and Laplace 
BVPs with the same boundary data. For exterior problems, it was shown in [H^ that 
the convex hull of singularities of Helmholtz and Laplace BVPs are the same, assuming 
the same boundary data. Therefore, the location of poles in a Laplace problem can serve 
as a guide for the location of sources in a Helmholtz problem (see 67[ for an example of 
this idea). In particular, locating the Helmholtz sources so that they enclose the convex 
hull of singularities of the Laplace problem guarantees that all of the singularities of 
the Helmholtz problem are enclosed. According to the conjectures in [25|, which have 
recently been proved [30j . this guarantees exponential convergence. 
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Appendix A. Iterated rational fitting 

The IRF method we use is based on a variant of [i^ presented in [i^ . As explained 
in Section we seek a rational function P/Q that approximates W, and we do so by 
solving weighted linear least-squares problems of the from (|2T]) . Using the quadrature 
rules of Section [3. 2 [ we discretize the continuous norm in (I2ip and write the minimization 
problem as, 

K N 



mm 

a,b 



bnzl - anzlW{zk, 



(A.l) 



in which the unknowns are the monomial coefficients of P{z) and Q(^), denoted b and 
a, respectively. To avoid the trivial solution, a normalization condition, to be specified 
below, must be imposed on one of the polynomials. 

Let us rewrite (jA.ip in matrix form, using the notation 



A = diag(Ai,A2,...,AK), 
Q = diag Q{zi), Q{z2), ■ . 



W 

Vkn 



diag 



Wiz,),Wiz2),. 
n^0.1,...,N, 



.,W{zk) 



(A.2) 



We obtain, 



A^Q 



V, WV 



(A.3) 



The conditioning of the Vandermonde matrix, V, can be very poor (sampling points on 
the real line, for example), or very good (sampling points on the unit circle, for example), 
and this can affect the conditioning of the minimization problem (|A.3p . To deal with this 
problem, we construct two orthonormal bases, one for the column space of A = A^Q^^V 
and one for the column space of B = AsQ-iWV. This can be done with the Arnoldi 
iteration fH, H^, by noting that the columns of A form a Krylov subspace. 



K = span[d, Zd, Z^d, . . . , Z^d], 



(A.4) 
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where Z = diag (21, ^2, ■ • ■ , zk) and d = diag(A2 Q i), and similarly for the columns of 
B. Once we have constructed the orthonormal bases, 



colsp(P) = colsp(A), PI'P = I (A.5) 
colsp(Q) = colsp(B), q}q, = I, (A.6) 



we rewrite (IA.3I) in the new basis, 



min ||[P,-Q]c||, 



(A.7) 



c 



and impose the normalization condition C2n = 1, which corresponds to Q{z) being monic 
in the new basis. The solution of the minimization problem is then obtained by solving 
the over-determined linear system 



in the least-squares sense. In (jA.8|) . we denote the column corresponding to the A^th- 
degree polynomial basis function by Qat, while the rest of the Q matrix is denoted by 
Qi:]v_i. The numerical stability of (|A.7p is generally much better than that of (lA.Sp . as 
the unknowns are coefficients of two sets of orthonormal columns. 

Once (jA.3|) is solved, the zeros of Q[z), which are the poles for the next iteration, can 
be obtained without transforming back to the monomial basis. They are the eigenvalues 
of a generalized companion matrix j46j . 



where e^r is the N column of the order- A identity matrix and H is the upper Hessenberg 
matrix obtained from the Arnoldi iteration used to generate the orthonormal basis for 
colsp(B). 

Appendix B. Determining whether a pole is inside D 

When a pole is close to the boundary, or the boundary shape is complicated, it may 
be nontrivial to find out whether a pole is inside D or not. By restricting ourselves to 
boundaries described by piecewise rational fimctions, a solution to this problem is readily 
obtained. Using Cauchy's theorem, a pole z' is inside D if and only if 



This integral cannot be calculated exactly using the rational Gauss-Chebyshev rules 
because the Chebyshev weight function is absent. However, it can be calculated analyti- 
cally by breaking it up into rational-function sections, writing the integrands as rational 
functions, and using a partial-fraction expansion. 



[P, -Qi:Af-i]c = Qat 



(A.8) 



(A.9) 




(B.l) 
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